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In this Letter, we use a general renormalization-group algorithm to implement Propp and Wilson's 
"coupling from the past" approach to complex physical systems. Our algorithm follows the evolution 
of the entire configuration space under the Markov chain Monte Carlo dynamics from parts of the 
configurations (patches) on increasing length scales, and it allows us to generate "exact samples" of 
the Boltzmann distribution, which are rigorously proven to be uncorrelated with the initial condition. 
We validate our approach in the two-dimensional Ising spin glass on lattices of size 64 x 64. 



The Markov chain Monte Carlo method [l[ has devel- 
oped into a universal computational approach in many 
disciplines of science and engineering, and it remains of 
great importance in the field of statistical physics where it 
originated more than 50 years ago. Indeed, many difficult 
calculations in high-dimensional spaces can be expressed 
(more or less formally) as the calculation of expectation 
values of an observable O 
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To compute the ensemble average on the left of Eq. (JTJ) , 
the Markov chain [on the right of Eq. (Q}] passes from one 
configuration x t at time t to the next one, Xt+i, most of- 
ten in a way respecting the detailed balance condition. 
An ergodic Markov chain, which can eventually reach 
any configuration x from any other x, has the property 
to visit configurations with probability oc e~ pE ^ in the 
limit of infinite simulation time r S i m — > 00, where the 
time average indeed coincides with the ensemble aver- 
age. In all Monte Carlo calculations, the convergence 
toward the stationary values of means and (connected) 
correlation functions is exponential with, for example, 
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The correlation time r corr provides a crucial scale, be- 
cause only Markov chain Monte Carlo simulations that 
have run for times much longer that r corr yield useful re- 
sults and are essentially free of systematic errors caused 
by the initial condition. 

Data analysis routines allow to reliably estimate the 
correlation time (and rclatedly, the error of a Monte 
Carlo calculation) whenever r corr <C T s - lm (sec, for exam- 
ple, for a discussion of the "bunching method"). On 
the other hand, it is extremely difficult to ascertain that 
a Monte Carlo simulation has indeed converged, that is, 
that the essential requirement for its validity, r s ; m 3> T corr 
(or at least r S j m > T corr ), is satisfied. This difficulty is 
very prominent, for example in the simulation of disor- 
dered systems, where the multidimensional space of con- 
figurations x is extremely rugged, and where the Monte 
Carlo dynamics is governed by many different time scales 



(the longest of which is the correlation time). It is often 
impossible to assure the validity of a Monte Carlo calcu- 
lation in disordered systems and other fields without re- 
sorting to the comparison with alternative methods, like 
exact analytic solutions, power expansions, careful finite- 
size scaling, etc. Bhatt and Young address this point 
in their classic paper Q when discussing their conver- 
gence test "In practice, though, one needs a criterion 
which determines whether any error made in being not 
quite in equilibrium is acceptably small (...). We know of 
no "rigorous" such criterion but we have found that the 
following procedure works well in practice." 

More than a decade ago, Propp and Wilson [f| realized 
a major conceptual breakthrough: they showed how to 
reformulate the Markov chain Monte Carlo algorithm so 
that it generates "exact samples" that have no correlation 
with the initial configuration, not even an exponentially 
small one, as in Eq. ([2]) . The generation of exact samples 
is of greatest interest in many real- world applications: 
it would solve the above-mentioned problem because the 
samples are exactly in equilibrium (they carry no memory 
of the initial state) and because the criterion is rigorous. 

Unfortunately, Propp and Wilson's procedure, termed 
"coupling from the past", has been notoriously difficult 
to apply to complicated physical systems, because it im- 
plies, as we will discuss later, the monitoring of the entire 
configuration space of a physical system under the Monte 
Carlo dynamics. In the ferromagnetic Ising model above 
the Curie temperature, the entire configuration space can 
be monitored, very elegantly, by making use of a partial 
order of spin configurations ([5j, see also 0,01). Huber 
has presented an interesting method that, unlike the 
partial ordering approach, can be made to work for dis- 
ordered or frustrated systems Q, but only at very high 
temperature. 

In the present Letter, we show how to apply exact sam- 
pling to more complex systems than have been treated 
before, namely the two-dimensional spin glass on large 
lattices at low temperature. We use a general, yet rig- 
orous, method which monitors all the configurations in 
the whole system, not directly (because there are far too 
many of them) but through "patches" of initially much 
smaller scale. During the calculation we gradually in- 
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crease the size of the patches, not unlike what is done in 
rcnormalization-group calculations. At the latest stages 
of the calculation, the patches are of the same size as 
the lattice, and the configurations on patches correspond 
to configurations on the entire system. Our approach to 
exact sampling works in a wide range of temperatures, 
and for quite large systems. We are able to judge its 
efficiency by comparison with a "naive" method which 
performs a standard Monte Carlo simulation for a repre- 
sentative fraction of all the configurations. 

The two-dimensional spin glass does not have a phase 
transition at temperature 1/(3 = T > 0, but it is already 
a quite complex system, due to the presence of disorder. 
During recent years, the two-dimensional s pin glass has 



the past approach to exact sampling relies on the generic 
property of Markov chains to couple during their evolu- 
tion 



been a test bed for new algorithms [3, LL0|, lU|, and has 
given rise to controversies concerning the specific heat 
capacity at low temperature [1, Oil [12] . In this model, the 
convergence time of Monte Carlo calculations is difficult 
to estimate. 

In the heat bath algorithm, the spin <Ji(t) (on site i = 
1,...,N) is updated using a uniform random number 
Ti(t) =ran[0, 1]: 
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where hi(t) = J^j Jij a j(t) is the local field. The square 
lattice is bipartite. This allows us to update one entire 
sublattice simultaneously, in one "sweep" , from the spins 
on the other sublattice, using a vector of random num- 
bers T{t) = {Ti(t), . . . , Yjy(t)}. Two subsequent sweeps 
update the whole lattice. In the remainder of the Let- 
ter, we explain our algorithm for a unique instance of 
the two-dimensional Ising spin glass, defined by a spe- 
cific choice of {Jij = ±1} for nearest neighbors i and 
j on the 64 x 64 square lattice with periodic boundary 
conditions. We note that, in principle, any configuration 
a is described through one sublattice, that is, on 32 x 64 
sites. 

In the "coupling from the past" approach, one consid- 
ers the simulation as running between an initial config- 
uration, at t = — oo, and the present configuration, at 
t = (see Fig. [IJ. The configuration <r(t = 0) is an 
exact sample, because it results from an infinitely long 
Monte Carlo calculation. It is evidently impossible to 
perform an infinitely long simulation, but we may pick it 
up at an intermediate time, t — to < 0, where, in prin- 
ciple, the Markov chain could be in any one of the 2 N 
configurations, and determine cr(t = 0), if the Markov 
chain "couples" between t = to and t = 0. This means 
in our context that under the dynamics of Eq. ([3]) all 
the 2 N possible initial configurations cr(t = to) yield the 
same configuration after a finite number r coup of sweeps. 
If the chain does not couple, we need to complement 
the sequence of random numbers {T(to), . . . , T(— 1)} by 
values corresponding to earlier times. The coupling from 
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FIG. 1: Schematic representation of coupling from the past: 
The spin configuration <r(t = 0) is completely decorrelated 
from the initial configuration at t — — oo. The configuration 
c(0) follows from the 2 N configurations at t — to if the chain 
couples between t = to and t = 0. 

A lower bound on the "coupling time" T coup is obtained 
by checking that several randomly chosen initial con- 
figurations have evolved toward the same state starting 
from time to. We implement this procedure in a naive 
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FIG. 2: Weights vs time (in sweeps) of the heat bath algo- 
rithm for an instance of the 64 x 64 Ising spin glass at /3 = 0.5 
(configurations that couple add weights). At t — to = 631, 
the 2 18 randomly chosen initial configurations have coupled 
toward a single configuration. Up to time t ~ to + 330, many 
configurations with very small weights subsist. 

algorithm which yield a rigorous value of r coup only if 
the initial configuration comprise the entire configura- 
tion space. We apply the naive algorithm for the present 
instance of the 64 x 64 Ising spin glass with random num- 
bers {T(i ), T(< + f ), . . . }, for 2 18 initial configurations 
randomly chosen among the 2 64x32 ~ 3.23 x f 616 config- 
urations of the entire configuration space. Two different 
configurations cr(t) and <x(i) may coalesce at time t + 1 
and will remain the same from then on (see also [3]), 
so that the number of configurations <x(i) decreases with 
time. We may define the "weight" of a configuration er(t) 




3 



as the fraction of the original configurations at time to 
that have evolved towards <r(t) (see Fig. [2]). In our exam- 
ple, the weight rapidly concentrates on a few configura- 
tions, and on a single one after t — to = 631 < r coup . All 
weights ^ 2~ 18 in Fig. [2] are expected to remain essen- 
tially unchanged for other choices for the <r(to) and even 
for an ideal simulation taking into account all the possible 
initial configurations. At a difference with our renormal- 
ization group approach, the naive algorithm cannot yield 
a rigorous upper bound for r coup and cannot prove that, 
indeed, r coup = 631. We note that the coupling time de- 
pends on the realization of the Markov chain, that is, all 
the vectors {T(t ), T(t + T coup - 1)}. 

We now calculate an upper bound for r coup (with given 
random numbers {T(to), . . . }). The practical generation 
of exact samples is straightforward (lil [TBI . 

As shown in Fig. [3l each spin configuration on the 
entire lattice can be broken up into configurations on 
patches k = 1, . . . , N, that is, pieces of the lattice. We 
define a set Sk(t) of "configurations on patch fc" and the 
product 

n(t) = Si(t) ®S 2 (t) <g> ••• 

<g> SN(t) I (pairwise compatible). (4) 

At the initial time to, the set Sk{to) contains all possible 
spin configurations on the patch k so that Q(to) contains 
a superset of the 2 N configurations on the entire lattice. 
At later times t, the Sk(t) contains a superset of all con- 
figurations er(t), restricted to patch k. As indicated in 
Eq. the sets are "pairwise compatible". Pairwise 
compatibility of Sk (t) and Si (t) means that any configu- 
ration on patch k must agree in the overlap region k n I 
with at least one configuration (as in Fig. [3]) of Si(t) on 
the neighboring patch I. Pairwise compatib ility is easy 
to enforce by a procedure we call "pruning " UJ. It con- 
sists in eliminating all spin configurations in Sk(t) that 
lack a compatible configuration in Si(t). However, pair- 
wise compatibility is far from sufficient to assure that all 
elements in £l(t) are valid spin configurations. We may 
"assemble" two patches k and I by constructing config- 
urations on the patch k U I from each pair of compat- 
ible configurations of neighboring sets Sk(t) and Si(t). 
In practice, we can assemble all the patches into config- 
urations on the entire lattice only if the sets Si(t) are 
sufficiently small. 

The Monte Carlo algorithm updates a spin on site i as 
a function of its nearest neighbors. It follows that a spin 
configuration in Sk{t) (the set on patch k, of size m x m, 
at time t) allows us to determine the spin configuration at 
time t+1 for all the spins on the "center" of the patch, all 
the sites that do not touch its boundary. In the example 
of Fig. [3j the patches are of size 4x4, and their centers 
of size 2x2. More generally, the center of an to x m 
patch is of size (m — 2) x (m — 2). 

Many different configurations of Sk(t) yield, after one 



center 
A i — 



center 





+\+ 




+ 






+ 




+ 


+ 


- + + 






+ + - 




+ 




+ 


+ 


— + 


+ 


+ 


- + + 






— + 





patch k 



spin configs 



patch / 



FIG. 3: A spin configuration cr = {<7i, . . . , ctjv}, and the set 
of configurations Sk and Si on patches k and I. The overlap 
region k n / consists here of nine sites. The configuration on 
the center of a patch, after one sweep, depends on the patch 
alone. 



sweep, identical center configurations. The different con- 
figurations on nine neighboring (to — 2) x (to — 2) cen- 
ters can be assembled to form, after pruning, the sets 
S k (t + 1), for k = 1, . . . , N. The set S k (t + 1) remains 
a superset of all configurations <r(t + 1), restricted to 
patch k. Because of the coupling property of the heat 
bath algorithm, the size of the sets Sk(t + 1) for small 
values of r = t — t , is generally smaller than the size of 
Sk(t). After many sweeps, the coupling property is off- 
set by the loss of information during the assembly of the 
centers into the to x to patches, and the size of the sets 
Sk(t) starts to fluctuate around a constant value. At this 
point it becomes convenient to assemble patches of size 
(to + 2) x (to + 2), and to apply the above procedure to 
these larger patches. In this renormalization procedure, 
the small length scales are effectively "integrated out" 
by the coupling property of the Monte Carlo algorithm, 
and one is able to construct the configurations cr(t) on 
ever increasing length scales. The pruning and assem- 
bly of patches has been implemented in the PERL pro- 
gramming language using hashing tables and referencing 
procedures. The programs also transmits pairwise com- 
patibilities from time t to t + 1. This is highly efficient 
during the assembly. 



262144 



4096 



check point 



X 




o 
X> 

S 
a 



200 300 500 550 631 
time (sweeps) t— 1 

FIG. 4: Coupling of the 2 64x32 configurations of the two- 
dimensional Ising spin glass on a 64 x 64 lattice at /3 = 0.5 
followed through patches of increasing size (compare with 
Fig.HJ. 
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We now apply the renormalization-group algorithm for 
determining a rigorous upper bound for r coup to the same 
instance of the 64 x 64 Ising spin glass that was already 
followed with the naive algorithm in Fig. [5] We follow 
the evolution of the heat bath algorithm through patches 
and centers of size 6x6 and 4x4, respectively, starting 
with all the 2 18 = 262144 possible configurations on each 
6x6 patch at the initial time t = ty. We increase the size 
of the patches after each 50 sweeps (see Fig. [4]). At later 
times, we are able to assemble all configurations on the 
64 x 64 lattice from the patches. At time t — to = 330, 
for example, the patches are of size 18 x 18, and the 
complete assembly yields 24 configurations cr(t) of the 
64 x 64 system (see the "check point" in Fig. |4}. As a 
consistency check of our algorithm, we verified that these 
24 configurations include the 4 configurations present in 
the naive algorithm at this point (see Fig. [2]). The algo- 
rithm couples at time t — to = 631, that is, it reaches one 
configuration per patch. During the whole procedure, 
no configuration was dropped, and therefore, necessarily, 
the complete assembly of these patches (of size 30 x 30) 
yields a single configuration er(t + 631). In this example, 
we find that the upper bound for r coup agrees with the 
lower bound from Fig. [2l 

To conclude our discussion of the instance of the 64 x 64 
Ising spin glass, we note that our renormalization ap- 
proach has coupled after the same number of sweeps as 
the straightforward Monte Carlo calculation of the naive 
algorithm in Fig. [2] In this sense, the description of 
spin configurations through the overcomplete set Q(t) of 
Eq. (d) is efficient, even though all the 2 32x64 - 10 616 
configurations that comprise the entire state space of the 
system have been monitored. However, although both 
simulations have converged in 631 sweeps of physical 
time, it should be clear that our renormalization algo- 
rithm remains very costly in CPU time (for details see 
flil]). This was the price to pay in order to exhibit an 
exact sample. 

We have tested our procedure for the Ising spin glass 
at lower temperatures. At (3 = 0.60, we checked on many 
samples of the that we can still couple the Ising spin 
glass without problems. The number of configurations 
per Tn x ttl patch levels off after awhile, and the renor- 
malization procedure of Fig. [4] becomes crucial. We have 
also successfully tested the procedure in the ferromag- 
netic Ising model below the Curie temperature. In ana- 
lyzing the behavior of our algorithm, the naive algorithm 
of Fig. [2] (which generated a sharp lower bound in our 
spin glass example) becomes a valuable tool. 

In conclusion, we have presented a renormalization 
group approach to exact sampling, and applied it to a 
large instance of the two-dimensional Ising spin glass. 
Our approach allows us to control a huge number of con- 



figurations by means of patches whose sizes increase dur- 
ing the simulation. 

Our method can be expected to work generally for 
models with local interactions but does not rely on spe- 
cial properties, as the partial ordering of configurations, 
that only hold for severely restricted classes of models. 
The fact that the model is two-dimensional only simpli- 
fies the assembly of local information about the configu- 
rations (the patches) but is not required to make it work. 
We have successfully implemented the approach for the 
three-dimensional Ising spin glass fill ] at temperatures 
lower than can be handled with a previous method 0, H[ . 

It would be most exciting if exact sampling could now 
be applied to even more complex problems, as to the 
three-dimensional spin glasses around the transition tem- 
perature, where the question of whether a Monte Carlo 
simulation has converged is especially difficult to answer. 
As mentioned throughout this Letter, the "coupling from 
the past" approach allows one to draw exact samples 
(that are proven to have converged) and it puts all Monte 
Carlo simulations for which it can be applied on an ex- 
tremely solid foundation. In this Letter we hope to have 
made a crucial step towards the application of exact sam- 
pling to physically challenging problems. 
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